#include "Integral.hpp"
#include <iostream>

__global__ void precomputeKernel( float *sinValues, float *cosValues,
				  float minX, float maxX, float minY, float maxY, float step )
{
  
  int idx = threadIdx.x + blockIdx.x*blockDim.x;
  
  float x = idx*step + minX;
  float y = idx*step + minY;
  
  sinValues[ idx ] = sin( x*(float)(M_PI) );
  cosValues[ idx ] = cos( y*(float)(M_PI) );
  
}

__global__ void integralKernel( float *values, float *sinValues, float *cosValues,
                                float minY, float maxY, float step )
{
  
 
  int workSizeY = (maxY-minY) / step;
  
  int blockX = blockIdx.x % 125;
  int blockY = blockIdx.x / 125;
  
  int idx = threadIdx.x + blockIdx.x*blockDim.x;
  
  int xIndex = (threadIdx.x%16);
  int yIndex = (threadIdx.x/16);
  
  __shared__ float sCosValues[ 16 ];
  __shared__ float sSinValues[ 16 ];
  
  if( threadIdx.x < 16 )
  {
    sSinValues[ threadIdx.x ] = sinValues[ blockX*16+xIndex ];
    sCosValues[ threadIdx.x ] = cosValues[ blockY*16+xIndex ];
  }
  __syncthreads();
  
  values[ idx ] = exp( sSinValues[ xIndex ] * sCosValues[ yIndex ] ) + 1;
  
}

double computeIntegral( double minX, double maxX, double minY, double maxY, double step )
{

  int workSizeX = (maxX-minX) / step;
  int workSizeY = (maxY-minY) / step;
  int workSize = (maxX-minX) * (maxY-minY) / step / step;
  
  float* devValues;
  float* devSins;
  float* devCoss;
  
  cudaMalloc( (void**)( &devValues ), workSize*sizeof( float ) );
  cudaMalloc( (void**)( &devSins ), workSizeX*sizeof( float ) );
  cudaMalloc( (void**)( &devCoss ), workSizeY*sizeof( float ) );
  
  precomputeKernel<<< (workSizeX+255)/256, 256 >>>( devSins, devCoss, minX, maxX, minY, maxY, step );
  cudaThreadSynchronize();
  integralKernel<<< workSize/256, 256 >>>( devValues, devSins, devCoss, minY, maxY, step );
  
  float* values = new float[ workSize ];
  
  cudaMemcpy( values, devValues, workSize*sizeof( float ), cudaMemcpyDeviceToHost );
  cudaFree( devValues );
  cudaFree( devSins );
  cudaFree( devCoss );
  
  double sum = 0.0;
  
  for( int i=0; i<workSize; i++ )
    sum += values[ i ];
  
  sum *= step * step;
  sum /= (maxX-minX) * ( maxY-minY );
  
  return sum;

}